#ifndef MESH_
#define MESH_

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cmath>
#include <time.h>

#include "DGMatrix.h"
#include "DGIntMatrix.h"
#include "ANN/ANN.h"
#include "Parser/Parser.h"

using namespace std;

enum elmtType {Point=15,Line=1,Triangle=2,Quadrangle=3};
//------------------------------------------------------------------------------
class Mesh{

private:
  
  int nbNodes;
  int nbGlobalElmts;
  int nbIntFaces;
  int nbBFaces;
  
  DGMatrix nodes;                                                               // (nbNodes,2)           --> nodes coord (x,y)
  DGMatrix Hmaxmin;                                                             // (2,2)                 --> [hminmin,hminmax;hmaxmin,hmaxmax]
  
  //by element type
  
  vector<int>               nbElmt;                                             // [EType]               --> nb elements for this EType
  vector<int>               nbElmtFaces;                                        // [EType]               --> nb faces for this EType
  vector<DGIntMatrix>       elmtNodesTag;                                       // [EType](eID,nbEnodes) --> nodes tag for eID
  vector<vector<int> >      elmtGmshTag;                                        // [EType][eID]          --> Gmsh Tag for eID
  vector<vector<DGMatrix> > elmtFaceID;                                         // [EType][eID](nbF,L&R) --> [eFID1,eFID2,... ; LorR1,...] for Tri=3, Quad=4
  vector<vector<DGVector> > elmtH;                                              // [EType][eID](2)       --> [hmin(inscribed),hmax(circum)]
  vector<vector<DGMatrix> > elmtNPcoord;                                        // [EType][eID](nbNP,2)  --> nodes coord (x,y) for eID
  vector<vector<DGMatrix> > elmtJ;                                              // [EType][eID](2,2)     --> eID jacobian 
  vector<vector<DGMatrix> > elmtInvJ;                                           // [EType][eID](2,2)     --> eID inverse jacobian
  vector<vector<double> >   elmtDetJ;                                           // [EType][eID]          --> eID determinant jacobian (detJ = area)
  
  
  //by interface
        
  vector<int>      bGmshGeoTag;                                                 // [bID]                 --> Gmsh geo Tag of boundary bID
  DGIntMatrix      bFacesNodesTag;                                              // (bID,2)               --> nodes tag for bID
  DGIntMatrix      bFacesTag;                                                   // (bID,4)               --> [ETypeL,EIDL,FNumL,BCTag]
  vector<DGVector> bFacesNormal;                                                // [bID](3)              --> [Nx,Ny,Norm]	(detJl = normal/2)
  DGIntMatrix intFacesTag;                                                      // (fID,6)               --> [ETypeL,EIDL,FNumL,ETypeR,EIDR,FNumR]
  vector<DGVector> intFacesNormal;                                              // [fID](3)              --> [Nx,Ny,Norm]
  
  //divers
  vector<DGIntMatrix> nodesIndex;                                               // [EType](fNum,node)    --> index of nodes for fNum of EType
  
private:
  
  void   buildJacobian   ();
  void   computeElmtSize ();
  
public:
    
  Mesh(ClassParameter&);
  ~Mesh();
  
  vector<elmtType> elmtGmshType;                                                // list of Gmsh elmt type
    
  int           getNbNodes         ()                            {return nbNodes;}
  int           getNbGlobalElmts   ()                            {return nbGlobalElmts;}
  int           getNbIntFaces      ()                            {return nbIntFaces;}
  int           getNbBFaces        ()                            {return nbBFaces;}
  
  int           getNbElmt          (int eType)                   {return nbElmt[eType];}
  vector<int>*  getNbElmt          ()                            {return &nbElmt;}
  int           getNbElmtFaces     (int eType)                   {return nbElmtFaces[eType];}
  vector<int>*  getNbElmtFaces     ()                            {return &nbElmtFaces;}
    
  int           getElmtNodesTag    (int eType,int eID,int nI)    {return elmtNodesTag[eType](eID,nI);}
  int           getGmshTag         (int eType,int eID)           {return elmtGmshTag[eType][eID];}
  
  double        getIntFacesNormal  (int fID,int i)               {return intFacesNormal[fID](i);}
  DGVector*     getIntFacesNormal  (int fID)                     {return &intFacesNormal[fID];}
  double        getBFacesNormal    (int fID,int i)               {return bFacesNormal[fID](i);}
  DGVector*     getBFacesNormal    (int fID)                     {return &bFacesNormal[fID];}
  
  int           getElmtFaceID      (int eType,int eID,int f)     {return elmtFaceID[eType][eID](f,0);}
  int           getElmtFaceID      (int eType,int eID,int f,int LR) {return elmtFaceID[eType][eID](f,LR);}
  int           getIntFacesTag     (int fID,int i)               {return intFacesTag(fID,i);}
  int           getBFacesTag       (int fID,int i)               {return bFacesTag(fID,i);}
  int           getBGmshTag        (int fID)                     {return bGmshGeoTag[fID];}
  
  int           getEGmshType       (int eType)                   {return elmtGmshType[eType];}
  double        getNode            (int nID,int dim)             {return nodes(nID,dim);}
  double        getHmaxmin         (int i,int j)                 {return Hmaxmin(i,j);}
  double        getElmtH           (int eType,int eID,int minmax){return elmtH[eType][eID](minmax);}

  DGMatrix*     getNPcoord         (int eType,int eID)           {return &elmtNPcoord[eType][eID];}
  DGMatrix*     getJ               (int eType,int eID)           {return &elmtJ[eType][eID];}
  DGMatrix*     getInvJ            (int eType,int eID)           {return &elmtInvJ[eType][eID];}
  double        getDetJ            (int eType,int eID)           {return elmtDetJ[eType][eID];}
  
  double dus           (int,int,int);
  double monomeof      (int,double,double,int);
  double dmonomeof     (int,double,double,int,int);
  
};
#endif
